Spontaneous plaquette dimerization in the J1 — J2 Heisenberg model 
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We investigate the non magnetic phase of the spin-half frustrated Heisenberg antiferromagnet on 
the square lattice using exact diagonalization (up to 36 sites) and quantum Monte Carlo techniques 
(up to 144 sites). The spin gap and the susceptibilities for the most important crystal symmetry 
breaking operators are computed. A genuine and somehow unexpected "plaquette RVB", with 
spontaneously broken translation symmetry and no broken rotation symmetry, comes out from our 
numerical simulations as the most plausible ground state for J2/J1 — 0.5. 
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The nature of the non magnetic phases of a quantum 
antiferromagnet is a topic of great interest and has been 
a subject of intense theoretical investigation since Ander- 
son's suggestion about the possible connections with 
the mechanism of high- Tc superconductivity. 

Within the Heisenberg model the simplest way in 
which the antiferromagnetism can be destabilized is by 
introducing a next-nearest-neighbor frustrating interac- 
tion leading to the so called J1—J2 Hamiltonian 
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where Si — {Sf , Sf , S^) are s — 1/2 operators on a square 
lattice. Ji and J2 are the (positive) antiferromagnetic su- 
perexchange couplings between nearest and next-nearest- 
neighbor pairs of spins respectively. In the following 
we will consider finite clusters of N sites with periodic 
boundary conditions (tilted by 45° only for N — 32) . 

Although there is a general consensus about the disap- 
pearance of the Neel order in the ground state (GS) of the 
present model for 0.38 < J2/J1 < 0.60 §-|, no definite 
conclusion has been drawn on the nature of the non mag- 
netic phase yet. In particular an open question is whether 
the GS of the J1 — J2 Heisenberg model is a resonating 
valence bond (RVB) spin liquid with no broken symme- 
tries, as it was originally suggested by Figueirido et al. 

The other possibility is a GS which is still SU{2) in- 
variant, but nonetheless breaks some crystal symmetries, 
dimerizing in some special pattern |^-|ll| . 

In this paper we address this point using exact diago- 
nalization (ED) and a quantum Monte Carlo technique, 
the Green function Monte Carlo (GFMC), which allows 
the calculation of GS expectation values on fairly large 
system sizes {L < 144). This is extremely important to 
draw reasonable conclusions on the physical thermody- 
namic, zero temperature, properties of the model. 

For frustrated spin systems as well as for fermionic 
models, quantum Monte Carlo methods are affected by 
the so called sign problem that can be controlled, at 
present, only at the price of introducing some kind of 
approximations, such as the fixed node (FN) one ||l^. In 



this work we have also extensively used a recently de- 
veloped technique, the Green function Monte Carlo with 
Stochastic Reconfiguration (GFMCSR), which improves 
systematically the accuracy of the FN approximation for 
GS calculations 

The FN method allows to work without any sign prob- 
lem by using the following simple strategy: the exact 
imaginary time propagator e"'^-^ - used to filter out the 
GS from the best variational guess j-i/'c) ~ is replaced by 
an approximate propagator e""^^"'" such that the nodes 
of the propagated state e~'^^™\ipG) do not change, due 
to an appropriate choice of the effective FN Hamiltonian 
Hfn (which in turn depends on IV'g))- The FN approxi- 
mation becomes exact if the so called guiding wavefunc- 
tion \iPg) is the exact GS. However for frustrated spin 
models even the best variational wavefunction of the Jas- 
trow type used to guide the FN dynamic, provides 
rather poor results even for the GS energy expectation 
value |4|jl|,|l|. 

The GFMCSR method allows to release the FN ap- 
proximation and to obtain results much less depend- 
ing on the quality of the guiding wavefunction. Dur- 
ing each short imaginary time evolution t ^ t + At, 
where both the exact and the approximate propagation 
can be performed without sign problem instabilities, the 
FN dynamic is systematically improved by requiring that 
a given number p of mixed averages |]l3| of correlation 
functions are propagated consistently with the exact dy- 
namic. By increasing the number of correlation functions 
one typically improves the accuracy of the calculation 
since the method becomes exact if all the independent 
correlation functions are included in the stochastic re- 
configuration (SR) scheme. 

Typically |4p^, few correlation functions {p ^ 10) 
allow to obtain rather accurate values of the GS en- 
ergy with an error much less than 1% - for lattice sizes 
{N < 36) where the exact solution is available numer- 
ically - and without a sizable loss of accuracy with in- 
creasing size. Such accuracy is usually enough to repro- 
duce some physical features that are not contained at the 
variational level, as it has been shown in a previous study 
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of the present model Q. In the latter work, in fact, with 
a gapless guiding wavefunction, it has been possible to 
detect a finite spin gap in the thermodynamic limit for 
J2/J1 > 0.4. 
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FIG. 1. Size scaling of the energy gap to the first 5 = 1 
spin excitation obtained with the GFMCSR technique for 
J2/J1 = 0.38 (full triangles), 0.45 (full squares) and 0.50 (fuU 
circles). Data for the unfrustrated (J2 = 0) Heisenberg model 
taken from Ref. [Q, are also shown for comparison (empty 
circles). Lines are weighted quadratic fits of the data. 

We have extended the previous GFMCSR calculation, 
with the same guiding wavefunction of Ref. |^ , by includ- 
ing in SR conditions not only the energy and all SfS^ in- 
dependent by symmetry, but also the antiferromagnetic 
order parameter. The latter, as discussed in Ref. p3| , 
though not improving the accuracy of the calculation, 
allows a very stable and reliable simulation for large p. 
The new results, extended up to iV = 144, confirm the 
previous findings of a finite spin gap for J2/J1 ^ 0.40 
(Fig. 11). 

As suggested in Refs. in order to investigate 

the possible occurrence of a spontaneously dimerized GS 
displaying some kind of crystalline order, we have cal- 
culated the response of the system to operators breaking 
the most important lattice symmetries. This can be done 
by adding to the Hamiltonian (|^) a term 50, where O is 
an operator that breaks some symmetry of H . On a finite 
size, the GS expectation value of O vanishes by symme- 
try for (5 = and the GS energy per site has corrections 
proportional to 5^ as by the Hellmann-Feynman theorem 
-de{S)/d5 = {d)s/N. Therefore e{S) ~ eo - X(5V2 , 
being x the generalized susceptibility associated to the 
operator 6, namely, x = 2{tpo\0{Eo - H)~'^6\tpo) /NJi. 
For ^ 00, if true long range order (LRO) exists in 
the thermodynamic GS, an infinitesimal field 6 ^ 1/N 
must give a finite {0)s/N ~ implying that the finite 
size susceptibility x = {0)s/5N has to diverge with the 
system size . Thus susceptibilities are a very sensitive 
tool for detecting the occurrence of LRO. 

We have considered the response of the system to the 
following symmetry breaking operators 



Oc = ^ (Si • Si+a- — Si ■ Si+y) , (2) 

i 

Op = ^e'Q°-^'S, -8,+, , (3) 

i 

with X = (1,0), y = (0,1) and Qo = (7r,0), for 
the rotation and the translation symmetry, respectively. 
Within ED and GFMC technique the susceptibility x — 
~d^e{S)/dS^\s=o can be evaluated by computing the GS 
energy per site in presence of the perturbation for few val- 
ues of S, and by estimating numerically the limit S ^ 
of the quantity xi^) = ^2(e(5) — eo)/S^. 
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FIG. 2. ED results for the J1 — J2 chain: Xpi^) associated 
to the operator Op (breaking the translational invariance) 
for J2/J1 — 0.2 (a) and J2/J1 — 0.4 (b). Data are shown for 
iV = 12, 14, 16, 20, 24, 30 for increasing values of Xp{S). 

As we have tested in the one dimensional Ji — J2 
model, the numerical study of LRO by means of x{^) 
is very effective and reliable. Here a quantum critical 
point at J2/J1 — 0.2412 separating a gapless spin fluid 
phase from a gapped dimerized GS, which is two- fold de- 
generate, is rather well accepted ]l9|-pT|. As shown in 
Fig. |, the response of the system to the perturbation 
50 p [Eq. (^], breaking the translation invariance with 
momentum fc = tt, is very different below and above the 
dimer-fluid transition point. However it is extremely im- 
portant to perform very accurate calculations at small 5 
to detect the divergence of the susceptibilities for large 
system sizes. 

In two dimensions, among the dimerized phases pro- 
posed in the literature, the so-called columnar and pla- 
quette RVB §-0 are the states which have obtained the 
most convincing numerical evidences. Both the colum- 
nar and plaquette RVB break the translation invariance 
but only the latter preserves the rotation symmetry. As 
also suggested in a recent paper by Singh et al. jll], the 
appearance of a columnar state can be tested by using as 
order parameter the operator Oc defined in Eq. (H). As 
shown in Fig. ||, the ED results for iV = 16 and = 36 
indicate that the susceptibility associated with this kind 
of symmetry breaking, xci decreases with the system 
size. Using the GFMCSR, described before, we have ex- 
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tended the calculation up to iV = 64. The GFMCSR 
calculations, which reproduce pretty well the ED data, 
rule out clearly the columnar dimer order. 

The above result is in disagreement with the conclu- 
sions of several series expansion studies [7| p^Jll[ . How- 
ever, as stated in Ref. the series for xc are very 
irregular and do not allow a meaningful extrapolation to 
the exact result. In our calculation instead, even the ED 
results for < 36, are already conclusive. 
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3. Exact and GFMCSR calculation of xc{5) associ- 
Oc (columnar dimerization) for J2/J1 = 0.5. 



This can be obtained by applying a generalized Lanc- 
zos operator (1 -|- aH) to the variational wavefunction 
\iPg)i where a is a variational parameter. This defines 
the so called one Lanczos step (LS) wavefunction, which 
has been particularly successful for the t- 
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-J model 

In the present model by using the LS wavefunction, 
a clear improvement on the variational estimate of the 
GS energy is obtained. More importantly, as shown in 
Fig. ^a), the LS wavefunction allows a much better es- 
timate of the susceptibility. Remarkably, on all the finite 
sizes where ED is possible, the GFMCSR estimate of 
this important quantity is basically exact within few er- 
ror bars (see also Fig. This calculation was obtained 
by including in the SR conditions the energy, the spin 
spin correlation functions up to next-nearest-neighbors, 
distinguishing also S'f 5'| and (5'f % + SfS^) (p = 4). 
The mixed averages of these correlation functions can be 
computed over both the wavefunction I'tpc) and the LS 
wavefunction {l + aH)\ipG) during the same Monte Garlo 
simulation. Thus with a LS wavefunction one can also 
easily double the number of constraints that are effec- 
tive to improve the accuracy of the method {p = 8). In 
this case we have tested that it is irrelevant to add fur- 
ther long range correlation functions in the SR conditions 
even for large size. 



Having established that the columnar susceptibility is 
bounded, it is now important to study the response of 
the Ji — J2 model to a small field coupled to the per- 
turbation Op of Eq. (^ , breaking the translation invari- 
ance of the Hamiltonian. The evaluation of xp, with a 
reasonable accuracy, is a much more difficult task. In 
fact in this case the ED values of the susceptibility for 
iV = 16 and = 32 increase with the size and much 
more effort is then required to distinguish if this behav- 
ior corresponds to a spontaneous symmetry breaking in 
the thermodynamic limit. As it is shown in Fig. ^(a), 
the FN technique, starting from a guiding wavefunction 
without dimer order, is not able to reproduce the ac- 
tual response of the system to Op, even on small sizes. 
The GFMCSR technique allows to get an estimate of the 
susceptibility which is a factor of three more accurate, 
but not satisfactory enough. In order to improve this 
estimate, we have attempted to include in the SR con- 
ditions many other, reasonably simple, correlation func- 
tions (such as the spin-spin correlation functions • Sj 
for |ri— Tjl > v^), but without obtaining a sizable change 
of the estimate of xp- In fact the most effective SR con- 
ditions are those obtained with operators more directly 
related to the Hamiltonian [p3 14|. 

After many unsuccessful attempts, we have realized 
that it is much simpler and straightforward to improve 
the accuracy of the guiding wavefunction itself. In fact it 
is reasonable to expect that both the FN and the GFM- 
CSR will perform more efhciently with a better l^pc), i-e., 
with an improved initial guess of the GS wavefunction. 
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FIG. 4. Xp{S) associated to Op (plaquette dimerization) 
for J2/J1 = 0.5, N = 32 (a), iV = 64 (b) and = 100 (c): 
FN (empty squares), GFMCSR (full squares), FN with LS 
(empty circles), GFMCSR with LS (full circles), exact (empty 
triangles) . 

By increasing the size, the response of the system is 
very strongly enhanced, in very close analogy to the one 
dimensional model in the dimerized phase (see Fig. ||(b)). 
This is obtained only with the GFMCSR technique, since 
as shown in Fig. ^ the combination of FN and Lanc- 
zos step alone, is not capable to detect these strongly 
enhanced correlations. For N = 100 the GFMCSR in- 
creases by more than one order of magnitude the response 
of the system to the dimerizing field. This effect is par- 
ticularly striking, considering that the starting guiding 
wavefunction is spin wave like |2^], i.e., gapless, Neel or- 
dered and without any dimer LRO. This suggests that 
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all our systematic approximations are able to remove al- 
most completely even a very strong bias present at the 
variational level. 

We believe that the numerical results we have pre- 
sented here give a very robust indication of a spontaneous 
dimerization with broken translation symmetry but with- 
out broken rotation symmetry (as discussed before), i.e., 
a plaquette RVB. This kind of state can be thought of a 
collection of rotation invariant valence bond states 




where |o o) — , Such plaquettes cover only 

one half of the possible elementary plaquettes of the lat- 
tice since two plaquettes cannot have a common side. In 
this way one necessarily has to break translation invari- 
ance and the resulting GS is fourfold degenerate in the 
thermodynamic limit, in agreement with the Haldane's 
hedgehog argument described in Ref. [p4| . 



—de{5)/dS ~ (e(0) — e{6))/S. As shown in the inset of 
Fig. H, the finite size effects of this quantity seem to satu- 
rate above the = 64 lattice size for 5 > 0.04, allowing a 
rather convincing estimate of the dimer order parameter 
as Op ^0.1, being sizably non zero. The sharp crossover 
of the size effects for > 64 is due to the presence of 
a finite triplet gap in the excitation spectrum (Fig. |^), 
implying, typically, a finite characteristic length. The 
value of the order parameter Op is rather large consider- 
ing that J2/J1 = 0.5 is very close to the transition point 
for the onset of sponaneous dimerization J2/J1 — 0.40. 
This is an interesting and measurable physical property 
that can be, in principle, investigated experimentally. 
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HTCS and LOTUS) and MURST (COFIN99). We thank 
A. E. Trumper, M. Calandra, M. Capone, F. Becca, G. 
Santoro, A. Parola and E. Tosatti, for fruitful discussions. 
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FIG. 5. Exact (empty triangles) and GFMCSR (circles) 
calculation of X-p(<5) (plaquette dimerization) for J2/J1 — 0.5 
and (from the bottom) A'' = 16, 32, 36, 64. Inset: numerical 
determination of the order parameter (see text). Lines are 
guides for the eye. 

In the past, among several attempts to guess the nature 
of the non magnetic phase of this model, the description 
closest to ours was that proposed by Zithomirski and 
Ueda . Amazingly, part of their conclusions were based 
on an unfortunate mistake in the series expansion pd| ]. 

The quantitative estimate of the order parameter can 
be obtained by taking first the thermodynamic limit 
A'' ^ 00 of the order parameter Op{S) = {0)s/N at 
fixed field 5, and then letting 5^0, lim Op(5) = Op 

being the value of the order parameter. In order to es- 
timate Op{5) at fixed size we have used the Hellmann- 
Feynmann theorem with a finite difference estimate of 
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